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ABSTRACT 


A set of spherical sampling functions is defined such that they are related to 
spherical -harmonic functions in the same way that the sampling functions of informa- 

2 

tion theory are related to sine and cosine functions. An orderly distribution of (N+ 1) 

2 

sampling points on a sphere is given, for which the (N+ 1) spherical sampling func- 
tions span the same linear manifold as do the spherical-harmonic functions through 
degree N. In this case, the transformations between the spherical sampling functions 
and the spherical -harmonic functions are given by recurrence relations. The spherical 
sampling functions of two arguments are extended to three arguments and to nonspherical 
reference surfaces. Typical applications of this formalism to geophysical topics are 
sketched. 


RESUME 

On dgfinit un ensemble de fonctions sphgriques telles qu'elles sont relives aux 
fonctions sphgriques-harmoniques de la m£me manilre que les fonctions d'echantillon- 
nage de la thgorie de 1 1 information le sont aux fonctions sinus et cosinus. On donne 

p 

uie distribution ordonn£e de (N + 1) points d'^chantillonnage sur une sphere, pour 

2 

lesquels les (N + 1) fonctions sph^riques d'echantillonnage recouvrent le m£me es- 
pace lin£aire que le font les fonctions sph^riques-harmoniques jusqu'au degre N. Dans 
ce cas, les transformations entre les fonctions spheriques d'echantillonnage et les 
fonctions spheriques-harmoniques sont donn^es par des relations de recurrence. Les 
fonctions spheriques d'echantillonnage pour deux arguments sont etendues cl trois ar- 
guments et cl des surfaces de reference non spheriques. On ebauche des applications 
typiques de ce formal isme aux problfemes geophysiques . 
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KOHCIIEKT 


OnpeaeJiHeTCH cwcTeua c$epwqecKHx npotfHbix cpyHKijHM TaKOBbix ^to 

OHM C Bfl 3HH bl CO CCjjepM yec KMMM- rapMOHHWeC KKMH (JyHKUMHMM T3KHM ie 
o6pa30M K.aK m npofiHbie (Jy hkhmm Teopnn MH$opMapHK c BH3 aHbi c $yHKLi,M- 

hwim CHHvca h KOCMHyca. npMBOAMTCH opraHH30 BaHH o e pacripeKeJieHwe 

2 2 

( N + l) npo^Hbix ToieK no c$epe, ajih KOTopbix (n + i ) c$epMHecKne 

npo^Hbie (JyHKHMM oxBaTM BaiOT to see canoe riMHewHoe MHoroo6pa3we khk 

m c$epwnecKMe rapMOHHTecKue $yHKU,ww CTeneHH . B stom cnyuae npe- 

o6pa30BaHne Mexny c$epwuec km.mm npo^HKMM c&yHKpnHMn m c$epMyecKMMM 

rapMOHH^eCKHMM $yHKl(MHMM npKBOflMTCH C nOMOIlIBK) peKypeHTHbIX OTHOUieHMM. 

C$epwyecKHe npofiHbie cjbyHKHMH HByx apryMeHTOB pacnpocTpHHflKTCH hjih 

Tpex apryivieHTOB m Ha Hec$epMHecKHe KOHTponBHbie njioma«M. HaPpacbiBa- 

KTCH TMnMHHbie npHMeHCHMH 3TOTO $OpMaJlM3Ma K. r eO$M3 K<iec KHM TenaM . 
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SAMPLING FUNCTIONS FOR GEOPHYSICS 


G.E.O. Giacaglia and C.A. Lundquist 


1. MOTIVATION AND HISTORY 

The mathematical formalism in the following sections arose initially as a candi- 
date technique for the analysis of satellite-to-ocean altitudes (Lundquist, Giacaglia, 
Hebb, and Mair, 1969). Many authors have noted the potential geophysical value of 
accurate measurements of the height of a spacecraft above the ocean (e.g., Frey, 
Harrington, and Von Arx, 1965; Godbey, 1965; Moore, 1965; Lundquist, 1967a), 
and various review panels have long recommended implementation of such measure- 
ments in the space program of the United States (e.g. , Pierson and panel, 1965; 
Stewart and panel, 1965; Woollard and panel, 1966; Doyle and panel, 1969; Kaula, 
1970) . Indeed, the program of the National Aeronautics and Space Administration now 
does include flights of spacecraft altimeters during the 1970s (Pierson and Mehr, 1972; 
Stanley, Roy, and Martin, 1972). Because preparations for analysis of the altitude 
data may require as much lead time as does fabrication of the spacecraft hardware, we 
were motivated to investigate appropriate analytical tools. 

The formalism summarized here evolved over several years. In 1967, we first 
examined the requirements for the analysis of satellite-to-ocean altitudes (Lundquist, 
1967a, b). This led to a judgment that an alternative mathematical representation of 
the geopotential would facilitate treatment of altitude data (Giacaglia and Lundquist, 
1968; Lundquist et al. , 1969). The latter reference defined in general terms a set 
of functions, differing from spherical harmonics, in which the earth's potential can 


This research was supported in part by grant NGR 09-015-002 from the National 
Aeronautics and Space Administration. Contract N000 14-67-0013 from the Office of 
Naval Research supported Giacaglia' s work at the University of Texas, Austin. 
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be expanded. The most important feature of this formalism is the property that the 
coefficients in the expansion are exactly a tabulation of the values of the potential at 
a grid of points on the earth. 

In the 1969 documents, several important mathematical topics were not yet 
treated. One such topic is the realization that the functions defined in those refer- 
ences are a natural extension of the sampling functions of information theory 
(Lundquist and Giacaglia, 1972a). Another related topic is the utilization of the 
sampling-function technique for geophysical quantities other than the gravitational 
potential (Giacaglia and Lundquist, 1971). 

The main objective of the present work is a comprehensive mathematical devel- 
opment of the sampling-function formalism for representing geophysical quantities 
around a spherical or a nearly spherical body such as the earth. The discussion 
here unifies and largely supersedes the mathematical developments in the earlier 
references. The final section concerns some geophysical applications of the 
sampling-function formalism that augment the applications discussed in the 1971 
references. 

Motivated by the satellite-altimetry expectations or by other circumstances in 
which a spherical-harmonic representation is awkward, various authors have offered 
still further alternative representations for the geopotential. These include distri- 
butions of point masses (Sjogren, Muller, and Gottlieb, 1971; Balmino, 1972; 

Kaplan and Kunciw, 1972; Needham, 1970) and surface layers of varying mass density 
(Vinti, 1971; Croopnick, 1970; Koch and Morrison, 1970; Koch and Witte, 1971; Koch, 
1972). For a critical review of alternative representations, see Kaula (1970). 

The principal motivation for an expansion in sampling functions instead of 
spherical harmonics arises when a representation containing much short -wavelength 
detail is required. If, for example, physical features of 1° extent are to be repre- 
sented in the geoid, a spherical -harmonic representation through degree and order 

180 must be employed. The value of the radius to the geoid is then calculated at any 

2 

point as the sum of (180+ 1) = 32,761 terms. If features of 0.°1 extent are considered, 

2 

the spherical -harmonic representation would have (1800+ 1) = 3,243,601 terms! With 
the equivalent sampling-function expansion, only a few terms contribute significantly 
to the value at any point. 
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It is a fundamental property of the sampling functions that they exist in one-to- 
one correspondence with a grid of sampling points on a sphere or similar reference 
figure. At each of these sampling points, only the single corresponding sampling 
function has a nonzero value, and this function multiplied by its coefficient makes 
the dominant contribution to the sum in the neighborhood of the sampling point. Thus, 
the coefficients of the terms in a sampling-function expansion are a tabulation of the 
value of the represented quantity at the grid of sampling points. 

On the other hand, the spherical-harmonic expansion through some degree and 

order N and the sampling-function expansion that spans the same linear manifold are 

2 

completely equivalent. In both cases, there are (N + 1) independent functions and 
terms in the truncated series. The spherical -harmonic functions and the sampling 
functions are related by a linear transformation in the linear manifold. Further, the 
transformation and its inverse can both be expressed analytically for any N. Hence, 
it is simple to transform to a spherical-harmonic representation when that is most 
convenient, or to a sampling-function representation when the tabular characteristic 
is useful. 
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2. GENERAL PROPERTIES OF SAMPLING FUNCTIONS 


2. 1 Introduction 


In this section, we address the concept of sampling functions appropriate to the 
representation of periodic functions over a well-defined domain. The main interest 
is in linear, spherical, and nearly spherical domains. 

To achieve a reasonably self-contained discussion and to establish notation, it is 
convenient to restate the definitions of well-known quantities and operations. (The 
reader is referred to the mainy standard references — for example, Sansone, 1959.) 

A function f(x) that is square integrable (in the Lebesque sense) over a certain 
domain Djx} will be indicated by 

f(x)GL 2 (D) . 

All the L 2 functions in D constitute a Hilbert space H(D). As usual, in this space 
we denote the scalar ( inner) product of f, g E L 2 P) by 

(f, g) = I f(x) g(x) dx , (2.1) 

J D 

which is a real number since we shall consider only real valued functions. Also, the 
positive number 

Ilf II = (f,f) l/2 (2.2) 

is the norm of f. The distance between two functions f, g £ L 2 (D) is the norm of f - g. 
In general, if the distance is small, the difference f - g need not be small everywhere 
in D. 
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A sequence of functions j f }£ L 2 (D), finite or countably infinite in number, is 
orthonormal in D if 


( f i> f k> = 6 ik 


for all i, k. The series 


(2.3) 


c 1 f 1 (x) + c 2 f 2 (x) + • • • + c n f n (x) + • • • 


(2.4) 


where 


c k ( f ’ f k> » 

is the generalized Fourier series of f(x). The c^ are the generalized Fourier coeffi- 
cients of f(x) with respect to the set fp f 2 , . . . . In general, nothing can be said about 
pointwise convergence (with respect to the defined norm) of series (2.4) to the function 
f(x), but it can be shown that (2.4) has the least-squares deviation (distance) from f(x) 
and that the Bessel's inequality 

£ cl * ||f | 2 (2.5) 

n 

holds. 

The sequence { f R J of orthonormal functions is complete in L 2 (D) if, for every 
f(x) £ L (D), the corresponding generalized Fourier series converges to f in the mean, 
i. e. , 


lim 

n -* ao 


n 




= 0 . 


( 2 . 6 ) 


6 



The Parseval relation 



k=l 



( 2 . 7 ) 


is a necessary and sufficient condition for completeness . Other conditions are also 
available, but this will suffice for our needs. It should be noted that if f, g (EL^D) 
possess the same Fourier series with respect to a complete system, then they are 
equal (almost everywhere) in D. 

Uniform convergence with respect to an orthonormal set implies uniform 
convergence (almost everywhere) to a function f. 

It will be necessary to consider the concept of completeness for nonorthogonal 
systems . If | f | £ L 2 (D) is a sequence of functions, linearly independent in D ( i. e. , 
every finite number of terms of the sequence is a linearly independent system), the 
sequence is complete in D if any f £L 2 (D) can be approximated in the mean by linear 
combinations of f R to any accuracy. 

Another concept of importance is the following: A subset n of a space S is called 

a linear manifold in S if it has the property that when the elements 1, , f n f_ £Q, 

— ~ ■ — ■ x a n 

any linear combination c-.f.. + c 0 f 0 + • • • + c f £ fi. IfQis closed in S, is called a 
J 1122 nn ’ 

subspace of S. If Q is a subspace of L 2 (D), then each element can be uniquely 
represented in the form 

f = g +h , 

where g Efl and where h is orthogonal to all elements of fi. The element g is called the 
projection of f into Q . 

If there is an m such that j f^, f 2 , . . . , f j is the maximum number of orthogonal 
functions, linearly independent in Q, then 0 is a finite dimensional linear manifold, 
or the span of { fp f 2 , . . . , f^J . 
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With respect to the problem of fitting a function with a series or a polynomial 
(which is actually a projection operation), it is important to remember the funda- 
mental theorem: Let L , f 0 , . . . , f , . . . E L 0 (D) be a given orthonormal sequence and 
let f EL 2 (D). The generalized Fourier series of f(x) is (2.4), where 


c 


k 



( 2 . 8 ) 


Let a positive integer N and an orthonormal system 


f 1 (x),f 2 (x),...,f N (x) 

be given. Then, from all functions of the form 


k l f l + k 2 f 2 + 


+ k N f N 


J 


the one that has the least distance from f is series (2.4), with the c^'s given by equa- 
tion (2. 8) . It is of paramount importance to note that this theorem holds, in general . 
only for orthonormal systems . 


For a finite dimensional linear space, the concept of completeness loses its 
meaning, in the sense that any maximal set of linearly independent functions is a basis 
of the linear space. In this case, the necessary and sufficient condition for the func- 
tions | fp f 2 , - • - , f n }■ to be linearly independent in D is the nonvanishing of the Gram's 
determinant 


G = !(f i ,f j )| , (2.9) 

but, of course, this condition also applies to linear independence of any subset of an 
infinite sequence j f n j . 

Finally, it is important to observe that, although the inner product has been de- 
fined without any weighting function, introducing a weighting function does not affect 
any of the previous statements. Weighting functions should nevertheless E L 2 (D) and 
be nonnegative in D. 
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2.2 Sampling-Function Concepts 


Before we deal with the concept of sampling functions, it is helpful to describe 
their origin, which can largely be identified with the theory of interpolation and 
Lagrange polynomials. From this point of view, the literature is quite abundant on 
one-dimensional or n-dimensional spaces R n , but very sparse on other types of spaces. 


Let S(D) be an N-dimensional subspace of a Hilbert space H(D), and fp f^, . . . ,f^ 
be an orthonormal basis in S. Any function f £ S is uniquely written as a linear com- 
bination of fp f^, . . . , f„. On the other hand, any g £ H is best fit in S by the projection 
of g on S(D) over the orthonormal basis j fp fg, ; . f^| , so that for any set of real 
numbers k x ,k 2 , ... ,1^, 




k. f. 
1 3 


( 2 . 10 ) 


on D (e.g. , Sansone, 1959; Rektorys, 1969). 

This is the classical problem of interpolation. Generalizations to a best fitting 
of f and its derivatives through a certain order are widely discussed in the literature, 
including some early papers by Hermite (1878). 

Now, consider another set of functions cpp cp^, . . . , cp^ £ S, such that 

/ 

N 

9j ( x ) = 5Z C ji f i (X > » j =1 , 2 >---> N , (2- 11 ) 

i=l 

and let C = j j be a non singular matrix. By definition, 



J 


(cp i jCp j) ^2 C ii {i V^ ^2 C i^ C ^j 


1=1 


1=1 


so that the matrix 

{(cp.,cp.)} = C T C (2.12) 

is nonsingular. It follows that cpp cpg, . . . , cp^ is a basis for S. 

Now let Q(x 1 , x 0 , . . . , x ) be a point of D, the domain of definition of all g £H(D). 
If the functions cp^(Q), k = 1, 2, . . . , N, can be defined in such a way that for a set of N 
points Q. £ D we can verify 

\(Qj) = 6 kj , (2- 13 ) 

then such functions are classically called collocation functions and are particular 
cases of sampling functions for the subspace S. In the case S = H, we shall have a 
denumerable infinite number of sampling functions. The completeness of these with 
respect to all g £ H(D) cannot be assumed a priori . 

The functions f^ £ S are uniquely represented everywhere in D by 


f k ^2 ’ k » 

j=l 


where {a kj } = {cj 


-1 


Since the a, . are real constants, 
kj ’ 


N 


f k (Q i> = £ a k J t Pj< Q £> = • 

j=i 


so that 


10 



N 

f k (Q) = X f k (Q j ) cp j (Q) 

pi 

uniquely. Consider any other function f £S. It can be written uniquely as 


k f k (Q) b kSWV Q) 

k=l k=l j=l 

N N N 

= I] Z {b k f k (Q j^ Vj (Q) = £ f(Q i } V P) ’ 

j=l k=l j=l 

so that, for all f £ S, Q £ D, 


N 

f(Q) = f (Qj) 'Pj(Q) , 
j=i 


(2.14) 


which is the basic justification for the name of sampling functions given to the cp.(Q). 
Analogously, any g £ H(D) has an expansion in S given by 

N 

g S (Q) = Y, g(Q j } ’ Q ED , 
j=l 


but since the set of cp. is not, in general, orthonormal in the integral sense, such an 
expansion cannot be assumed to be the best fitting of f in the sense already defined. 
It is obvious that, in an N-dimensional subspace, the N conditions (eq. (2. 13)) com- 
pletely define cp^. Therefore, orthogonality of these, in the sense 


C? j ) 


cpjjx) cp.(x) u>(x) dx = 0 , 


(2.15) 
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for k ± j, can be achieved, if possible at all, only by a particular choice of the sampl- 
ing points Qj (j = 1, 2, . . . , N) G D or a proper choice of a weighting function oj(x) >0 
in D. In some instances, orthogonality can be obtained or will follow automatically 
as a basic property of the subspace. 

In general, sampling functions satisfy the following properties: 

A. Orthogonality over the set of sampling points, 


N 


w vv 


= 6 . 


ij 


k=l 


(2.16) 


B. Normality over the set of sampling points, 

£ ‘W = E ^ m 1 ■ <2 - 17> 

k=l j=l 

C. Normality everywhere in D, i. e. , 

N 

^(Q) = 1 , Q GD . (2.18) 



This last property follows immediately from the fundamental property (2. 14). In fact, 
suppose 

N 

y] \(Q) = f (Q) G S(D) ; 
k=l 


then 
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N 

f(Q) = ^ f (Qj) cp.(Q) , 

and, therefore, f(Q^) - 1 for j = 1, 2, . . . , N, so that, necessarily, f(Q) = constant = 1. 


2. 3 Sampling Functions and Fourier Analysis 

It is not necessary to discuss in this report the nature and properties of sam- 
pling functions and how they are defined in their most general form. However, in 
view of the problem we originally proposed to solve, it is instructive to mention 
the role of these functions in information and communication theory. Historically, 
sampling functions are related to a generalization of the process of interpolation 
introduced by E. T. Whittaker (1915) under the name of cardinal functions. The 
basic statement, sometimes referred to as the Shannon principle but actually proved 
by J. M. Whittaker and several others, is the following: If a function f(t) contains 

no frequencies higher than q; radians per second, f(t) can be recovered from periodic 
samples fjt^) taken 2rr/2o; seconds apart . More general results are available, such 
as the recovery of f(t) from nonperiodic samples, deletion of a finite set of samples, 
conditions for obtaining f(t) when only the past is sampled, and a minimum mean- 
square error estimate of f(t) based on any (eventually nonperiodic) set of samples 
(Beutler, 1961). The function f(t) can be a random process, a periodic process 
affected by noise, or a function generally supposed to be in a Hilbert space . 

Functions containing frequencies not higher than a certain value are called band- 
limited functions. 


The nature of our original motivation is contained in the above description — 
namely, the recovery of the mean sea level from samples corresponding to satellite- 
borne altimeter measurements. Nevertheless, we are facing a much more general 
problem, for no sampling functions, in the above sense, have been introduced for 
functions defined on a sphere, spheroid, or ellipsoid. Take, for example, the mean 
sea level; it can be represented as a series of spherical harmonics: 
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(2.19) 


r <C,X) = R 0 £ £ + b nm V nm (C,X)] , 

n=0 m=0 

where the X nm and Y nm are defined in equation (2. 52). If r(£, X) E Lg (sphere), the 
series (2. 19) will converge in the mean. Moreover, in truncating the series to some 
degree n = N, i. e. , writing 


N n ~ \ 
r ( £ X) = R n V V a X +b Y 
N' ’ ' 0 2-~t ' 11111 nm nm nm / 

n=0 m=0 


( 2 . 20 ) 


we know that this projection on the linear manifold defined by the finite set of spheri- 
cal harmonics of maximum degree N has the property of least-squares deviation from 

r(£, X) if the coefficients a and b are defined by 
’ ' nm nm J 


a nm 4 nK 0 J_i< r( ^’ X) X nm (£ > X) dX 
and 

1 f 1 f 2rr 

b = A P dC| r(C,X)Y (£,X)dX . 

nm 4ttR q •'-i *'o nm' ’ ' 


( 2 . 21 ) 


( 2 . 22 ) 


In a generalized sense, r^(£, X) is a band-limited function in two variables, 

C = cos 9 and X, which are, respectively, the geocentric colatitude and longitude. In 
more detail, let a grid of longitude and latitude circles be drawn on a sphere. To 
this grid, there correspond a minimum and a maximum value of wavelength both in 
longitude and in latitude. 


In a geoid representation, the grid spacing should be chosen to correspond to the 
resolution of geoidal undulations attainable with the observation instrument. This is 
pertinent in the case of satellite-borne altimeters or in the case of the resolution of 
harmonics above the noise level in the perturbations introduced by the geopotential in 
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the orbits of satellites. These measurements correspond, respectively, to 0? 1 wave- 
length for the geoidal undulations ultimately measured by an altimeter and about 10° 
wavelength for the geopotential variations sensed by satellite orbits. 

In the usual terminology of information theory, the quantities we want to describe 
— mainly the geoid and the geopotential — are treated as band-limited functions 
both in longtiude and in latitude and therefore can possibly be recovered from samples 
taken at maximum intervals corresponding in some sense to half the wavelength. 
Specific relations corresponding to these quantities will be given in Section 6, where 
applications are discussed. 

In one dimension, the mathematical correspondent to the Shannon principle can be 
stated simply. If f(t) is band-limited, i.e. , if 

/ +Q 

e“ luJt g(u>) duo , (2. 23) 

-0 


where g(uu) is integrable, then f(t) has the series expansion 


= ^ f sin(flt - nTT ) 
1(t) 2^ l nt - nn 

n=-°o 


(2.24) 


On the other hand, if f(t) is periodic, it is evident that equation (2.24) reduces to a fi- 
nite number of terms. If f(t) is not periodic, the sampling theorem can be used from 
the point of view of interpolation; i.e., we require that, given a suitably restricted 
sequence of numbers, we construct a function f(t) defined by the (finite) series (2.24) 
that interpolates values between these numbers. Evidently, in doing so, we reverse 
the situation and force f(t) to be a band-limited function, for its spectrum will vanish 
outside the frequency interval [— Q , +fi] . Interpolation in these terms and correspond- 
ing bounds for truncations in the one-dimensional case are well known (e.g. , Beutler, 
1961; Zakai, 1965). We only point out here that bounds can be defined for an 
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interpolatory process based on sampling functions, and in some instances the sequence 
of truncated series 


V 

can be proved to converge (at least in the mean) to f(t) as N— 00 . 

The connection between the sampling theorem and Fourier analysis is easily 
established. Consider the infinite series 


+N 

■E 

n=-N 


f ( nn N sin (Qt - nrr) 
\ Q / Qt - nrr 


f(t) (a cosnuot+b 

' ' 2 0 ' n n 

n=l 


sin nu>t) , 


(2.25) 


which converges under the Dirichlet conditions on f(t). Let f^(t) denote the partial 
sum 




1 N 

= 7 T a„ + V' (a cos ntut + b sin mot) , 
2 0 2-j n n 


n=T 


which can be written 


r rr/u) 

f(T)K^[»(T-t)]dT . 

-TT/ (JO 


where 


KnOO 


sin f(2N + l)(x/2)l 
2 sin (x/2) 


(2. 26) 


(2.27) 
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If we extend the interval 2-rr to the entire real axis, it is readily verified that equa- 
tion (2.27) reduces to sampling functions, and equation (2.26) to the corresponding 
theorem. For periodic functions, they represent the same result. 

Finally, it is important to observe that the most general form of the sampling 
theorem corresponds to random processes that are band-limited, in terms of ex- 
pected values, while the form of the theorem remains the same with the condition 
that f(t) is a distribution. Assuming the response to a pulse to be affected by a ran- 
dom noise, the application of the sampling theorem should also provide a method of 
recovering the pulse by assuming band limits and therefore eliminating noise 
problems (Campbell, 1968; Zakai, 1965; Balakrishnan, 1957; Bennett, 1956). 


2. 4 Classical Examples of Sampling Functions 


A typical example, and probably one of the oldest, is given by Lagrange interpola- 
tion polynomials. Consider all polynomials R(x) of maximum degree N defined on the 
segment a < x < b. Let x. ( i = 0, 1, . . . , N) be N + 1 points on the segment, with 
x Q = a and x^. = b. The pertinent sampling functions are Lagrange polynomials: 



" N 


N 

V x) = 

TT < x - x i ) 

/ 

TT <x k - x i> 


i=l 

/ 

i=l 


i^k 

/ 

i^k 


k = 0, 1, . . . , N . (2.28) 


On the other hand, L k (x) are not an orthogonal set, and in fact, as is well known, they 
do not satisfy the condition for best fitting of L£(a, b) functions. 


This idea can be generalized by considering all trigonometric polynomials T(X) 
of maximum degree N, defined on the line segment 0 ^ X < 2n. If X. ( i = 0, 1, . . . , 2N) 
are 2N + 1 points on the segment, the pertinent sampling functions are the 2N + 1 
trigonometric polynomials 
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p k (M = 


“ 2N 

j [* sin 2 (X " 

/ 

r 2 b 

n 

' sin | (X k - X.) 

i=0 

/ 

i=0 

_ i/k 

/ 

i/k 


In the particular case where the X. are uniformly distributed, 


X. = 


2rri 


"i 2N + 1 ’ * 0,1,..., 2N , 

we obtain the well-known sampling functions 


(2.29) 


k = 0, 1, .. . ,2N 


(2. 30) 


\ (X) 2N + 1 


N 


1+2 cos j (X - X^) 
1=1 


1 sin{[(2N+ l)/2](\ -X k )} 


2N + 1 sin (1/2) (X - X^) 


, k = 0, 1, . ,.,2N . (2.31) 


All these functions and their properties are described in detail in the literature (see, 
e. g. , Brillouin, 1962, pp. 93-111). They can easily be generalized to two or higher 
dimensions, for squares, cubes, or hypercubes, both for polynomials of several 
variables and for multiperiodic functions. It is interesting to note that the q k (\) are 
also orthogonal in the integral sense: 

( \' q j )= / 0 'itfiTr'kj ■ 


The functions 


\< x > 


(2.32) 
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are a set of orthonormal sampling functions in [0, 2tt] and, therefore, are best fitting 
for any function f(\) £ L 2 (0, 2rr); that is, they satisfy an equality of type (2. 10). 

The examples we have presented so far are well established, and there is no need 
for specific derivation. Nevertheless, in general, given a linear manifold defined by 
the space of a finite number of functions and given an equal number of sampling points 
on the support of the manifold, there seems to be no broadly recognized definition of 
an equivalent basis constituted by sampling functions. This will be done in the following 
section. 


2, 5 Recursive Algorithm Defining a Finite Set of Sampling Functions 

Consider a set of functions f^Q), f 2 (Q), • • . , (Q) spanning an N-dimensional 

subspace L of a Hilbert space H(D), Q C D. Let a set of N sampling points 
Qf, . . - , £ D be given. Consider the problem of defining N functions 

S^(Q), S 2 (Q), . . . , S n (Q) spanning the same linear manifold and having the properties 
of being orthonormal over the set of points Q^, Q 2 , . . . , Q N , i. e. , 


N 


[ Si , Sj ] = £ S.(Q k ) S.(Q k ) = 5.. 
k=l 


(2. 33) 


Following the philosophy of the Gram-Schmidt process, let 


S, = 


“l 1 ! 

2. 

1„ 

o' L 

+ Ck'r-S, 

2 2 

2 1 

3. 


a„f„ 


3 3 

3 2 

k. 

k-1 

a, f, 
k k 

+ "k 

N, 

N- 

& N f N 

+ q, n 


\ s i 


4 s i 
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By imposing, in succession, conditions (2.33), we find 


j k 
oft = -o', 




'k " “k [sTTsT] ’ 

and therefore 


k = 1, 2, . . . , N; j = 1, 2, , 


(2.34) 


s k = “ k f k - 


k(. i£i[f k >s] 




k = 1, 2, . . . , N . 


(2. 35) 


The procedure is applicable, of course, if none of the inner products fS ^ , S^ ] is 

zero, which sets a restriction on the distribution of the sampling points Q , Q „ , . . . , Q .^ 
in D . If the distribution is nonsingular, can be defined so that [ S^, SjJ = 1, and 
therefore 



k = 1,2, . ..,N . 


k r t 

The are also obtained by recurrence by using [ S^, S k J =1, i. e. 


(2. 36) 


( a k) 2 j[ f k’ f kl -yj [ f k> S jJ 2 | = 1 ’ k = l,2,...,N . (2.37) 

For example, 

("l)Vi’ f i] = i 

(" 2 ) {f f 2’ f 2 ] " [ f 2> S ll } =1 
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The above method is quite important in the problem of constructing for a linear 
manifold a base that is orthogonal over a set of points, i. e. , orthogonal in a discrete 
sense. If orthogonality both in the integral and in the discrete sense can be achieved 
by a proper choice of a grid of points, we have an optimum solution to the problem. 

2 3 

Let us take, as a first example, the sequence 1, x, x , x , . . . , which is com- 
plete in Lg(a, b) for every closed interval [ a,b]. Orthogonalization leads to Legendre 
polynomials Pq(x), P^(x), P£(x), . • . , and the normalization can be set at will. Now 
consider the following discrete formulation of the problem, where the initial sequence 
is truncated at the Nth degree. In other words, we look at the linear manifold 
defined by 

2 N 

1, x, x , . . . , x , (2. 38) 

so that any polynomial of the Nth degree at most can evidently be written uniquely as a 
linear combination of the set (2. 38), or of the set Pq(x), P^(x), . . . , P^(x). Take N + 1 
equally spaced points x^ (j = 0, 1, . . . , N), with x Q = a and x^ = b, so that 

x. = a + ih 

1 

h = (2.39) 

j = 0,1, ...,N . 

One can construct the N + 1 polynomials ^q( 1,N), $^(j,N), . . . , § N (j,N) orthogonal over 
the set of points j = 0, 1< . . . , N, with the property that 


N 

[ * s (j,N),F s _ 1 (j)] 3 £ yj> N ) F s _i(3) = 0 > (2. 40) 

where [ $, F] now plays the role of an inner product, and F (j) is any polynomial of 
degree s - 1 or less. This leads to 
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= 


= 


2i_i 

N 

12.1 2 + 2N(N - 1) - 12Nj 
2N(N - 1) 


under the normalization 


N 


(j,N) $ /( j, N) = 6 / . 

n VJ> / n VJJ / nn 

3=0 


(2.41) 


The substitution 


j = ? (! + y) > 


which makes j = 0, (x = a), correspond to y = -1 and j = N, (x = b), to y = +1, leads to 


= 


= 


= 


p 0 (y, N > 

PjCy.N) 

pJ(y,N) 


= i 


= y 


3Ny - (N + 2) 
2(N - 1) 


(2.42) 


It follows that any polynomial of maximum degree N of y defined in [ -1, +1] can be 

* * * * 

uniquely written as a linear combination of P Q , Pp P 2 , . . . , P^. The important 
property is that for N — «>, it can be verified that 

P*(y, N) - P n (y) ; 

i. e. , the usual Legendre polynomials are obtained. This is an example of the con- 
struction of a base for a linear manifold from a discrete point of view, leading, in the 
limit, to a complete sequence in a Hilbert space. 
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As a second example, consider the sequence 1, cos x, sin x, cos 2x, sin 2x, . . . , 
which is orthogonal and complete in i, (-tt,tt). Any f(x) £ ,n) as a convergent 

(in the mean) Fourier series, as is well known. Moreover, if f(x) and f '(x) are 
piecewise continuous in [— rr , tt] , the series coincides with f(x) everywhere except at 
the points of discontinuity, where it gives the value (l/2)[f(x + 0) + f(x - 0)]. This 
produces the well-known Gibbs phenomenon. 

Again, consider a finite set 

1, cos x, sin x, cos 2x, sin 2x, . . . , cos Nx, sin Nx (2. 43) 


of 2N + 1 functions. In this case, the formulation for a discrete number of points is 
immediate since the set (2. 43) is also orthogonal over the set of 2N + 1 points, 


x. 

3 


• 2rr 
j 2N + 1 


j =0,1,...,2N 




and any trigonometric polynomial of maximum degree N, defined in [— tt , tt] , is 
exactly and uniquely given as a linear combination of the set (2.43). 


A third example, very important for our future problems, employs the sequence 


X 00’ X 10’ X ll’ Y ll’ X 20 J X 21’ Y 21’ X 22’ Y 22’ ' ' * 


5 


where 


X 

nm 

Y 

nm 


= P (C) cos m\ , 
nm' w ’ 

= P (C) sin m\ , 
nm' w ’ 


m = 0,l,...,n; n = 0, 1, 2, . . . 
m = 1 , 2, . . . , n; n = 0, 1 , 2, . . . , 


(2.44) 


are the usual surface spherical harmonics. They constitute an orthogonal and com- 
plete set in [-1< C — 1; 0< \ < 2tt] ; i.e., the unit sphere (a) for all functions 
f(C, \) £ 1 ^( 0 ). Any f of this space has a convergent (in the mean) Fourier series. 

In general, the discrete formulation of the present case is to our knowledge an un- 
solved problem. In fact, considering the linear manifold defined by the set 
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X 1 , X 0 , .... X „ of spherical surface harmonics of maximum degree n = N, we are 
1 1 (N+l) 2 2 

faced with the question of how to distribute (N+l) points on the sphere so that the 

Gram-Schmidt discrete orthogonalization process does not lead to zero divisors. 

More explicitly, assume that a point distribution Q.(£., X.), j = 1, 2, . . . , M (with 
2 3 3 3 

M = (N + 1) ), is given. As in the previous examples, discrete analogies of the above 

finite set spanning the same linear manifold can be defined by the relations (2. 33) 

through (2.37). In this specific case, 


S 1 = «1 X 1 
S 2 ~ a 2 X 2 + a 2 S l 


*k 


k._ , k-l„ . , 1 0 

= "k x k + 0, k s k-i + "' + “k s i 


(2. 45) 


o . _ M v , M-l 0 , , 1 c 

S M “M X M ®M S M-1 M S 1 ’ 


where, imposing the discrete orthogonality conditions over the set of Q^, 


I s 2 ' S l3 = 0 

[S 3 ,S 2 ] = [Sg.Sj = 0 


we find 


<*i 


-Of, 


[ x k» s t] 

rv^T ’ 


k = 1, 2, . . . , M; j = 1, 2, . . . , k-1 


(2. 46) 


Evidently, the set Q, is acceptable if the inner products [ S., S.] are nonzero for 
2 K .33 

j = 1, 2, . . . , (N+l) . The discussion of such problems is far from trivial, and we 

shall say more about them in Sections 3 and 4. 
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After the functions S^(Q) are obtained, the definition of the sampling functions is 
immediate. In fact, we look for a set of N sampling functions Wp Wg, . . . , 
spanning the same linear manifold and satisfying the conditions 

w j(Qk) = 6 jk • (2-47) 

If the W's span the same space as do the S's (and the f's), then, certainly, we can 
write uniquely 


N . 

W fe (Q) = PjJs.(Q) . (2.48) 

Multiplying equation (2.48) by S.(Q) and summing over all sampling points, we get 

[W k ,s.] = . 

3=1 

Since 

[W k ,Si] - S.(Q k ) 

and 

[Sj.Si] 

it follows that 

4 = S j( Q k ) , 

so that 

N 

W k (Q) = 2 S J (Qk) S j (Q) ’ 

3=1 

which solves the problem. Reciprocally, 


(2. 49) 


(2. 50) 
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(2.51) 


N 

S k (Q) = Yj S k (Q j } W j (Q) • 

1=1 


2. 6 Properties of Sampling Functions on a Sphere 

Since our objective is to define sampling functions spanning the same linear mani- 

o 

fold defined by the set of (N + 1) spherical harmonics up to a given degree N, we shall 
deal at present only with properties pertinent to this manifold. The problem of point 
distribution will be treated in the next section, so here we will describe properties 
derivable from the sole assumption that a set of sampling points exists, in the sense 
that the recursive algorithm described in Section 2. 5 is not singular. 


We shall eventually use spherical (surface) harmonics normalized according to 
the following relations: 


X (£, X) = 2 £ ' P (Q cos mX , m = 0, 1 . . . , n 

nm' s ’ ' V m nm' w ’ ’ ’ ’ 

(2.52) 

Y nm (C,X) = 2 P nm (0 sin mX , m = l,2, ...,n 

(see eq. (6. 6) for the relation to other normalizations) for n = 0, 1, 2, . . . , N, where 
£q = 1 and £ m = 2, m > 0 . Also, by the usual definitions, 


P (0 = 

nm' w 


2n + 1 (n - m) 1 
2 (n + m) ! 


1/2 


P (0 

nm' w 


and 




(i - c 2 ) m /2 d D+m 


2 n n! 


dC 


n+m 


(c 2 - 0" 


(2. 53) 


(2. 54) 


To simplify notation, let any of the spherical harmonics be represented by X, (£, X), 

2 K 
k = 1, 2, . . . , (N+l) . The necessary formulas to implement such single-index 
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notation will be given in Section 4. 1. The orthonormality conditions are written 


< X k’ X j> 



X k (C,X) X.(£,X) dX = 6 kj . 


For unnormalized spherical harmonics X , Y (or X, ), 

nm’ nm v k /5 

and 


we have (X , Y //) = 0 
' nm’ n m ' 


(X 

(Y 


nm’ 

nm’ 


X / 
n m 


Y / 
n m 



£ m (n +m)! 
(2n + l)(n - m) ! 


6/6 / 
nn mm 


(2.55) 


Any function contained in the linear manifold defined by the X k (£, X) can, by 
definition, be written uniquely as 



where 


a k = ifJ 0 d ^J 1 f(C ’ X) x k ^’ X)dc 

= (£,X k ) • (2.57) 

The definition of a pertinent set of sampling functions corresponding to the finite 

2 

number of wavelengths represented by the (N + 1) spherical harmonics rests on the 

2 

hypothesis that there exists at least one set of (N + 1) points (£ k , \ k ) such that (N + 1) 
functions W k (C, X) can be defined that satisfy 



and span the same linear manifold; 
should not vanish. 


that is, the Gram determinant G = [ (W, , W.) 

K J 


We must initially verify the relations 


_ M _ 

x k (C ’ X) = 53 w i (C,X) * (2,58) 

j=i 

2 

where M = (N + 1) and where the matrix 

should not be singular, so that equation (2. 58) can be inverted. It follows that 


where 


M _ _ 

W.^X) =£A jkXk (C,X) , 


k=l 


(2. 59) 


a = { A jk } ■ {\Y l • 

From equation (2.58), assuming X Q = P Q0 = 1, it follows that the normalization con- 
dition already mentioned holds everywhere on the sphere; i. e. , 


M 

^W.(C,X)=1 . (2.60) 

From equation (2. 59), it also follows that 

M _ 
k=l 
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and therefore 


j(W.,W.)j = AA T , (2.61) 

so that the Gram's determinant is not zero. Reciprocally, it follows from equation 
(2. 61) that if G ^ 0, | A I ^ 0 and therefore A ^ exists. 

When we use the unnormalized spherical harmonics, we shall define 

(X k ,X.) = c k 6 kj , (2.62) 

the inner product always being considered with the weight l/4rr, the area of the unit 
sphere . 

We shall now obtain a general form of the sampling functions W k (G, \) appropriate 

to the discussion of properties unrelated to the point distribution. We shall also freely 

alternate from the notation X, or (X , Y ) for unnormalized and X, or (X , Y ) 

k nin urn. k rim nm 

for normalized spherical harmonics, whichever happens to be more appropriate. At 
all times, M will mean (N + 1)^. 

From 

M 

W.(G,X) = 2 V X k (C ’ X) ’ 
k^l 

it follows that 

M 

< 2 - 63 > 

• i 

and therefore 

M 

W.(C,X) =£±-(W y X k ) X k (C,\) . (2.64) 

k^l k 
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Also, from 


M 

X k (C,X) = W.(C,X) 

it can be seen that 

M 

(Xr-'V = £ x k«.e,x 1 )<w 1> w J ) : 

i=i 

then, substitution into equation (2. 64) yields 

M M n 

W i <C ’ X) = E E 3T x k< c e’ X Jt )V j* X k (C > X) > < 2 - 65 > 

k=l (_X=1 k 

where 

= * 

Introducing the explicit expressions for X k (£, X) and cr^ into equation (2.65), we obtain 

N n £ (2n + 1) (n - m) ! M 

w i< c ’ k > ■ Z Z — — E V 

n=0 m=0 1 = 1 

x P nm<9 P nm (C) 003 m < X ‘ V • < 2 ' 66 > 

Let the cosine of the angle between the sampling point Q^(C^, ^) and a generic point 
Q(£, X) on the unit sphere be represented by y^, so that 

C V* " £ 2 ' cos " X ) • 
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By the addition theorem of Legendre, it follows from equation (2. 66) that 


M N 

W .(C,X) = £ v Jje £ (2n + 1) P n (Y,) 


1=1 n=0 


(2.67) 


or, in view of a well-known identity, 


1L 


M v. 

w i (CX) = (N + 1) E~ 

i-i 


P N^V “ P M^ (Y,) 


N+l' Y l , -I 


( 2 . 68 ) 


which is nonsingular at y^ = 1, since this is a simple root of the quantity in brackets. 

For computational purposes, equation (2. 68) is of no use since the quantities v . . are 

j Jo 

not known unless all are defined. Nevertheless, the equation is important for the 
study of the condition under which the point distribution leads to orthogonality with 
unit weight, that is, under what conditions 


(W j5 w k ) = ^ 6 kj 


6i.i . 


(2.69) 


If this relation were satisfied, equation (2. 68) would give 


w.(G, X) 


V j 1 - y. [ P N*V P N+l (V j ) ] ’ 


(2. 70) 


and after 




is imposed, it is found that 


P N ( V = P N+1 ( V ’ j^ = 1 >2,3,...,M 


(2.71) 
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The set of equations (2. 71) gives the necessary angular spacing, arccos y^, between 
the sampling points if the W's are orthogonal (in the integral sense) with unit weight. 
However, it does not give the location of the points on the sphere. The location is 
obviously not arbitrary, for we require the matrix to be nonsingular. Be- 

cause of the trivial solution y. . - 1, equation (2. 71) is actually an equation of the Nth 

}!L 

degree in y. „ . Since it can be shown that all roots are real, we conclude that a set of 
2 ^ 

(N + 1) sampling points leading to orthogonal sampling functions with unit weight can- 
not have more than N different angles between the points. 

2 

The simplest case, N = 1, corresponding to (N + 1) = 4 points, leads to a very 
interesting result. For N = 1, equation (2.71) gives 


_ 3 2 1 

2 y ]H~ 2 ’ 


(2. 72) 


whose solutions 'are y.^ = 1 (trivial) and y.^ = -1/3. The angle between any two of the 
four points must be the same and equal to about 109? 5. This is precisely a regular 
tetrahedron. Indeed, a tetrahedron with a vertex at the north pole leads to a set of 
four ort hogonal samplin g funct ions spanning the same linear manifold as 

(l, C. Vl - Q 2 ' cos X, Vl" - Z 2 ' sin X }. 


For N = 2, that is, nine points, the possible angles are two, whose cosines are 
given by y. . = (-1 + \/lP)/5. We find no configuration that satisfies such conditions. 
For N s: 3, the only practical way of solving the problem is by numerical means, but 
experiments in this sense have given no clue to a possible solution. Nothing can be 
asserted for general values of N. Moreover, it may be possible that by the introduc- 
tion of a weight, orthogonality can be achieved among the W's for a given point 
distribution. Another condition of orthogonality is easily obtained from the foregoing 
equations. Now assuming normalized spherical harmonics (o^ = 1), we find 


M 

£ = 


i=i 



(2. 73) 
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but this condition is of no use, because the v^'s are not known. It gives, nevertheless, 
a geometric interpretation to the constants : The v ^, are weights necessary to render 
the set of spherical harmonics orthogonal over the set of sampling points . 
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3. SAMPLING-POINT DISTRIBUTIONS 


3.1 General Considerations of the Point Distribution 


The problem of distributing in the most regular way a given number of points on 
a sphere has been studied for centuries and, in general, is still unsolved. A good 
account is given by Coxeter (1962), who discusses finding N(cp), the greatest 
possible number of small circles of angular radius c p that can be packed on the surface 
of a sphere without overlapping. Such a number has an upper bound, attained by only a 
few specific values of cp and given by 

N(cp) < N U(9) =^^76) » (3- D 

where sin a = (i/2) sec cp. A lower bound has also been established: 


o__ 

2/3 2 \ 3 2 

N(cp) > N t (cp) = 

cosec cp - L, cosec cp) 

L vT 

V / 


where N L is positive only for cp < 4l°48' (i. e. , cosec > 3/2). The most uniform 
point distribution is seemingly provided by the centers of the aforementioned circles. 
For N = 3, 4, 6, and 12, it is found that N = N^, with cp given by 60°, 54°44' (regular 
tetrahedron), 45°, and 31°43', respectively. 

In the problem of sampling functions for a sphere, the number of points has to be 

2 

a perfect square, (N + 1) , where N is the maximum degree of the spherical har- 
monics considered. This fact sets one important constraint. A second, and far 
more important, constraint is that the matrix |x^| = |x^(Cj, ?u)j should not be 
singular. 
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It is all too easy to construct point distributions for which jx^j is singular. 
Numerical experiments indicate that it is difficult to find an orderly, equatorially 
symmetric distribution that yields a nonsingular jx^j • No useful general charac- 
terization of singular or nonsingular cases emerged from our brief examination of 
this problem. However, the distribution discussed below seems to be one useful 
example of a convenient distribution of sampling points on a sphere. 


3. 2 The Recommended 1-3-5 Point Distribution 


o 

Lacking a general theory of how best to distribute (N + 1) points uniformly on a 

sphere, we wish to produce a distribution that is, in some sense, a generalization of 

the tetrahedron configuration. In one orientation of this configuration, one point is at 

the north pole and three points are on a circle south of the equator. Since spherical 

harmonics of maximum degree N are expected, by wavelength arguments, to correspond 

to a net of at most N + 1 latitude circles (parallels) and 2N + 1 longitude circles 

(meridians), it is natural to consider the pole and N different latitudes corresponding, 

2 

respectively, to £q, £ 2 , . . . , C N * A total number of (N + 1) points is obtained if we 

set one point at £q (north pole), three at five at C 2 , • • • > and 2N + 1 at £ N - This 
distribution will be called the 1-3-5 distribution. 

Now it must be decided what the values of Cp C 2 > • • • > Cvj are and h° w the 2k + 1 
points are to be distributed on each circle of latitude specified by £ Again, because 
of the natural wavelength in longitude \ of the spherical harmonics, apparently the 
points on each circle should be equally spaced so that at each £^ we associate 2k + 1 
points of longitudes, 


\j = i 


2rr 

2k + 1 




2k 


where ( ) indicates the integer part of the enclosed argument. This choice is made 

so that at subsequent latitudes we shift the longitude origin to 180°, and thus there is 
no accumulation of points at X = 0. In view of the results of Section 2. 6, regardless 
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of how the circles are situated and the points distributed on the circles, this 1-3-5 
distribution cannot with unit weight yield orthogonality in the integral sense among the 
W's. This can be illustrated by the case N = 2 discussed in Section 2, for which the 
cosine of the angle between sampling points must be ( -1 + \/(P)/5 for orthogonality. 

If 0^ and 0£ are the colatitudes corresponding to 

COS 9 1, 2 “ ” 5 T— 5~ = { 0.29 ’ 0^13377, 9 2 ~ 73. 1 , 

the cosine of the angular distance of three points on a circle at 6^ is given by 

2 2n 2 

cos <y = cos 0 ^ + cos -g- sin 0 ^ = 0. 214 , 

and that of five points on a circle at 0 2 by 

2 2tt 2 

cos 0 = cos 0 2 + cos — sin 02 = 0.367 . 

Thus, approximately, a = 78? 8, 0 = 68?5, which is significantly different from 0 2 -. 

If there is a general proof for the existence or nonexistence of a set of sampling 
points leading to orthogonal sampling functions, it would probably have to be 
based on an approach similar to that indicated in Section 2. 

The point distribution we propose here is given by 


_ , ,.k krr 

C k (-1) cos 2n + i > 


k = 0,1, ...,N 


- 0 - 1 


(3.3) 


2k 


where ( ) is the integer part of the argument, so that any point has coordinates 

( £ k , ^ k j). For this reason, we initially attach two indices to specify the sampling 
functions ; i . e . , 
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W k)< C ’ X) = E £[^r X nm 


— nm— 




n=0 m=0 


,»] 


(3.4) 


The distribution (3.3) is reasonably uniform, as shown in Figure 1. The cosine 

of the angular distance between points (£, , X, .) and (£ , X ) is given by 

K K] p pq 

/ 2 \ 1/2 / 2\ 1/2 

“Vpi'l'-fi) v c v) + ^ cos % - v • <3 - 5 > 

Figure 2 shows how the maximum and the minimum separations of two sampling 

o 

points change with (N + 1) and compares the upper and lower bound curves defined 
by equations (3. 1) and (3.2). 

Now there remains the problem of whether the 1-3-5 point distribution is every- 
where dense in the limit N — °°. Given a positive and arbitrary e and a point (G, X) op 
the sphere, consider a spherical rectangle with its center (0, and sides e. Then, 
for N sufficiently large, we should be able to prove that however small e may be, 
there is at least one point of the distribution lying inside the above defined rectangle. 


Consider < e/2 and e 2 < e/2. For some N, k, j, with 0 < k < N and 0 < j < 2k, 
we must seek to verify the following conditions: 


0 - e, < 0, <0 + 6, 

lk 1 

1 " e 2 < X kj < 1 + e 2 ’ 


(3.6) 


where 


0 ~ k 2TT 
k 2N + 1 ’ 


2rr 


k = 0, 1, . . . ,N 


(3.7) 


X kj ^ 2k + 1 + ^Nk 11 » J 0, 1, . . . , 2k 
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DEGREE OF REPRESENTATION 
WHERE TOTAL NUMBER OF POINTS = (DEGREE + I ) 


POINT SEPARATION 


Figure 2. Curve 1 is the maximum and curve 3 is the minimum spacing between 

points for the 1-3-5 distribution on a sphere. Curve 2 is the upper bound 
and curve 4 the lower bound on the number of points as a function of 
separation, as given by equations (3. 1) and (3.2), respectively. 



(For convenience in this proof, we use a slightly different subscript convention from 
that in equation (3. 3)). The coefficient is 0 or 1. It follows that k, j have to satisfy 
the inequalities 


9 - e- 9 + e- 

(2N + 1) < k < (2N + 1) 2 tt 


(3.8) 


and 


2 k + 1 
2 rr 


(X- 


P Nk n 


" e 2> 


<j<^(X + 9. 


Nk 17 + e 2> 


(3.9) 


under the restrictions (3. 7) on k and j. Given 0 and 9,X, the above conditions can be 
satisfied if N is sufficiently large and k, 3 are conveniently chosen. Actually, 9 + 0 ^ 
cannot be larger than rr, and ^ - 0 ^ cannot be smaller than zero, so that whatever 9,0 
are, we find 

0 < k < N +~ , 

which is always verified. On the other hand, through the same reasoning, 

0 < j < 2N + 1 , 

which is also within the allowed limits. 


Obviously, since k < N and j < 2k, we obtain from the relations (3. 8 ) and (3. 9) the 
limits 


(2N + 1) < N 

^ + P Nk + e 2 

( 2 k + 1 ) p - < 2 k . 


(3.10) 


The first of these gives, excluding the value 9 + 0 ^ = tt, 
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9 + e. 


2N > 


rr - (8 + 

and the second gives 

X + tt + e 


2N > 2k > 


rr - (X + e 0 ) 


so that 


2N > max 


X + rr + e. 


9 + e, 


rr - (X + e 9 ) n - (0 + e 1 ) 


(3.11) 


(3. 12) 


(3.13) 


Once N is found from equation (3.13), k can be chosen to satisfy relation (3.8), and 
thus condition (3.9) defines j. It is obvious that if 0 = rr and -* 0, then N -*«», as 
should be expected since the south pole is reached only in that limit. 

The distribution will be uniform in the limiting sense if the following is verified. 

Let n be the number of points falling inside an area AS of the sphere (unit radius), 
however chosen, and M be the total number of points. Then we must have 

lim n _AS 

.. M 4rr ' 

N -* co 

For the 1-3-5 distribution, the above uniformity condition is an open question. 
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4. SAMPLING FUNCTIONS FOR A SPHERE 


4. 1 Explicit Form of the Sampling Functions for a Sphere 

In Section 2 we gave the properties of sampling functions that are independent of 
the choice of a point distribution. Numerical tests indicate that the 1-3-5 distribution 
of Section 3 is nonsingular for any N tested, of which N = 36 was the highest. A gen- 
eral proof that it is nonsingular for any N is at present not available. Nevertheless, 
the nature of the 1-3-5 distribution makes it unlikely that it would lead to a singular 
matrix mainly because of its rather complete lack of symmetry. Of 

course, this is a weak mathematical justification. With this reservation in mind, we 
will proceed to derive expressions for the sampling functions of two arguments cor- 
responding to the 1-3-5 distribution of sampling points on a sphere. 


In order to obtain the explicit form of the functions W, .(£, X), written as 


N n 

v c - x,= e e 


n-0 m=0 


A" m X (C,X)+b!™Y (C,X) I , 

kj nm' ’ ' kj nm' s ’ 7 J ’ 


— nm — 

’kj nm' 


(4.1) 


it is convenient to proceed in the following steps: 

A. Reduce all double-index coefficients to a single -index notation. 

B. Solve the recurrence relations defining the intermediary set of "spherical 
harmonics" S p (£,X) that satisfy the condition [S p ,Sq] = 6 pq (as in Section 2.5). 

C. Obtain single-index coefficients corresponding to equation (4. 1) from these 
recurrence relations and from the definition of the W's. 

D. Revert to double-index coefficients and recover the A's and B's. 
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Step A. Single-coefficient specification of spherical harmonics and points of 
distribution. 

We first transform the set j X nm , Y nm J to the set { X p | by the following relations 


where 


X - X , p=n+ 2m + 1 ; m = 0, 1, 2, . . . , n 

nm P ’ ’ ’ ’ 

— — 2 

Y _ X , p=n +2m ; m = 1, 2, . ..,n , 

nm p ’ ^ ’ ’ ’ ’ 


n = 0, 1, 2, . . . , N ; p = 1, 2, . . . , M ; M = (N+1)‘ 


Conversely, 


(4.2) 


X 


nm 


P 1 y 


nm 


where 


<a/F> - P 


^ 0 , Case 1 below 


= 0 , Case 2 below. 

Again, ( ) indicates the integer part of the argument. 


(4.3) 


For Case 1, 


n = (n/p" 1 ) ; 


1 2 

even: m = (p - n ) — Y 


p - n = 


2 

odd: m = i ( p - n 2 - 1) — X. 


nm 


nm 


For Case 2, 


n -,/iT - 1— X 

nn 


Next we transform the set X^) into (^ , X q ) by the relations 
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q = k 2 + j + 1 


0 ( ^ 0 , Case 1 ' below 

<V7’> 2 -q{ 

1 = 0, Case 2 ' below , 

where for Case 1 ' 

k = ( > , j = q - k 2 - 1 , 

and for Case 2 ' 

k = /cp -1, j=q-k 2 -l. 

Then the single-coefficient sampling functions can be written 
M _ _ 

V c ’ X)= Z A qpV c ’ X) ' 

p=i 

We now have 

W q«s>V*‘q. >' s > q = 1 > 2 M • 

For notational convenience, it is useful to define the following matrices: 

W = { W qs} S K <C «’ X s>} = 1 (Me,ltlty) 

x ={v = v c s-v} 

A ={a 1 

1 qpj 

AX = I . 

Step B . Intermediary spherical harmonics. 

As we saw in Section 2.5, it is convenient to introduce the functions 

S P (C ’ X)= Z a pqV C ’ X) ’ 
q=l 


(4. 4) 


(4.5) 


(4.6) 


(4. 7) 
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with 


M 


[S,S 1=VS(C,X)S(C,X) = 6 
L P Q J 2s P s’ s' q Vto s’ s' pq 

s=l 


The recurrence relations defining the S's were obtained as 






j =1,2, .. .,p-l 





<* 3 S. 
P J 




a j S. 

P J 


- 1/2 


As an example, we list the first few of these spherical harmonics: 

S l = l 
S 2 = V^ 

S 3 = ~ £ 2 ' sin x 

s 4 = Z 1 " ? cos X 

S 5 = y/5~' ( 3 ^ 2 - l)/2 

sinX 

S ? = \/TP Cf/l - C 2 ' cos X 
S g = y/TUT (1 - e 2 ) sin 2X 
S 9 = y/TOT 1 (1 - Q 2 ) cos 2X . 


(4.8) 
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Step C . Single-index coefficients . 


All constants orj (p = 1, 2, . . . , M ; j = 1, 2, . . . , p) and S^(C k , X fc ) = S- k 
{j = 1, 2, . . . , M; k = 1, 2, . . . , M) are needed to implement the computation of the sampling 
functions W (£, X). These constants can be computed sequentially, as follows: 

r 


a\ - (N + l) -1 

s iq ~ (n + i) -1 , q = 1) 2, . . 


1 -1 M - 

*2 = “< N + 1 > L X 2s 

s=l 

“2 = )E(^s + “'2 S l s ) 2 j 

S 2 q = ‘4( 5 2q + “2 S 1,) 


M 


For i = 3, 4, . . . , M ; j = 1, 2, . . . , i-1 , 


and 


a? = 


M 


I =-E x - s - 


is js 


s=l 

M I _ i-1 . 

1 + 
s=l \ j=l 

i-1 . 


1 - 1/2 


S,_ = otU X,„ + ^ 


1/ x 

iq i| iq i jQ 

j=l 


In summary, the a's and S's can be generated sequentially in the following order: 
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1 

Of* 

1 

S 1P S 12’**’’ S 1M 
1 2 
01 2' 01 2 

S 21> S 22* ‘ * * » S 2M 
1 2 3 

’ a 3 * a 3 

S 31’ S 32’ ,,,J S 3M 

« 

* 

* 

12 k 
a k ’ a k ’ * * * ’ 

\v S k2’ ' * ‘ ’ S kM 
for k- 1,2, , ,.,M, 

By successive substitution of the recurrence relation 


V J=1 

into itself, we arrive at 

A n-1 n-k / n-k+1 n-k+2 

v«VE E E . E 

I k_1 J k 1 Vk-1 V'' 1 i k-2° J k +2 

** Hji ) 2 j 

- E "" >i “ii j 

JiV 1 - 1 

If we accumulate terms j^-1, j^=2, . . j^-n-l, it follows that the coefficients of 
S n (C,X) in terms of X^(C, X), as given in equation (4, 7), are 



(4.9) 
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(4. 10) 


n 

a = at 
nn n 


n-p / n-k+1 n-k+2 

a =«" ^^1 v y 

np n\ n p 2-i Zj 

k=2 Vk-l =P+1 J k-2 =pt2 


"li Jl >1 >2 

• • • \ at a. a. 

Lu n h h 

j^p+k-1 

These have the recurrence relation 


*? * p 
] k-l P 


n 

a = a V a . cP. . 
np PP L nj j 

j=P+l 


(4. 11) 


Step D . The sampling functions . 
As we have already shown, 

M 

q=l 

so that if we use equation (4.7), 


M ( M 

EE 

s— 1 ' q=s 


W P (C ’ X)= E E s qpVr s (C ’ X) ’ 


and therefore the coefficients A of equation (4. 4) are given by 
_ M 

' A qp = S s jq a jp 5 q= ;p = i,2, ..., M 

j=p 


(4. 12) 
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The index q corresponds to the pair (k, j ) , and p to the pair (n, m), so that 

Af m , B? m are easily recovered from A . 
kj ’ k] J qp 


4. 2 Structure of Sampling Functions for a Sphere 

The previous section gave explicit formulas for W, .(£, X) in a form appropriate 

KJ 

for computations. Another derivation of these relations follows, but in a form that 
exposes more of the general properties of the 

From equation (4. 1), it follows that 


n ( r n 


w k) (C k’ X) = E] e 5 kr p nm <v cosmx 

m=0 / n=m 


(C ) sin mX 
kj nm' k 7 


N / 

= e(*S 

—A ' 


cos mX + sin mX j . 


(4.13) 


We require that 


W kj^k’ X k^ 6 j 1 ’ 

which defines uniquely a set of 2k + 1 sampling functions (X) on the segment [0, 2tt] ; 
that is, according to equation (2.31), 


'V c k- x > = V’ 1 ’ = arn: 1 + 2 £ cos m <*• - V ■ 


k = 1, 2, . . . , N ; 
j = 0, 1, . . . , 2k . 

(4.14) 
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It follows that 


N 


V < p „o<V = 2irrr 

n=0 


N 


TnmTr 


E A, . P (Ci ) = oi ~ ~ T T cos m ^i • 
k] . nm'x 2k + 1 kj 


(4.15) 


n=m 

N 


E B kr p nm( £ k> = 2irri sinmX kj 


n=m 


for k = 1, 2, . . . , N ; j = 0, 1, . . . , 2k ; m = 1, 2, . . . , k ; and that 


N _ 

L a , 11111 p (C, ) = o 

kj nm'x 

n=m 


N _ _ 

E bT p (C ) = o 

k] nm'T 

n=m 


(4. 16) 


for k = 1, 2, . . . , N ; j = 0, 1, . . . , 2k ; m = k+1, k+2, . . . , N. 

Since C, k + 0, is not in general a zero of P (£, ), it follows that 
k nxxi K 


-nm =nm 
A, . = B, . 
kj kj 


= 0 


for m = k+1, k+2, . . . , N ; n = m, m+1, . . . , N. 


(4. 17) 


It follows further that we can write 


w tj (C,M = 


k j- 

z[ 

m=0 


T^(C) cos m\ + S™(C) sin mA. 


]■ 


(4.18) 
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where 


N 


T S (C) = t: A kj M P nm^^ 


n=m 

N 


(4.19) 


S Si (C >'E B k” P nm<S> • 


n=m 


(Note: S^(0 is not the same function as S p (C, X) of the previous section.) 

We next consider two situations leading to the definition and specification of addi- 
tional types of sampling functions appropriate for certain linear manifolds. 

We consider first the set 


l c c 2 c 3 C N 


or the equivalent set 


P 00 ,Pio(0,...,P M n(C) • 


NO' 


We wish to define the pertinent set of sampling functions corresponding to the points 


C k = (-l) k COS , k = 0, 1, . . . , N , 


(4. 20) 


that is, the values of Q defining the 1-3-5 point distribution. 


The sampling functions are easily found to be 


C k ( ^ 2N + 1 


1+2 E yv y<> 

3=1 


(4.21) 


where t.(£) = cos j 9 , Q = cos 9, and, as usual, = 1, = 2, k > 


0. 


52 



Consider next the set of associated Legendre functions of the first kind and of given 
order m: 


P 

mm 


(Q,P 


m+1, m 




N, m 


(C) • 


As is well known, this is an orthonormal set and defines a finite dimensional (dimen- 
sion = N - m + 1) space. Again, for the point distribution (4.20), we wish to define 
the pertinent set of sampling functions. The solution of this problem, which is not 
trivial in detail but is simple in essence, is given by application of the method 
presented in Section 2. 5. 


We begin by constructing a system orthogonal over the set of points C^, 
k = 0, 1, . . . , N. Evidently, there are m points in excess, which must be excluded; for 
reasons given below, we exclude the set £q, . . . , Specifically, the problem 

is to find, for any given N and m, a set of N - m + 1 functions, 


<w« = x;*r p „ m«> > (4 - 22) 

n=m 


that satisfy the N - m + 1 conditions 

G ^m ( 9 = 6 li > Z ’ 1 = m > m+1 > • • • » N • (4. 23) 

The intermediary functions with discrete orthogonality properties over the set 
C m’ C m+l'"-’ C N given by 

"mm '=■’ a mm 

%(0 

where a . is chosen so that 


-E[ P to>V] H )m<4 

j=m ' 


H = m+l, m+2, . . . , N , 
(4.24) 
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We can now write 


N 


G W£> - E h r H nm<0 > 


n=m 


multiply by H. m (Q, 


N 


V (C) G /m<0 - E h 7 H jm<° H nm<0 * 


n=m 


set £ = £ , m < s < N, and sum over s. This gives 
s 


N 


[V’ G «m] = E h " m ( V> H nml 


n=m 


or 


V (C *> - h f • 

and finally, 

G W» = E H nm«V H nm<» - < 4 ‘ 25 > 

n=m 


which, combined with equation (4.24), gives the desired result. In this way, we find 
the functions H nm (Q and therefore the G^ m (Q in a recurrent way. 


The form (4. 22) can also be obtained directly. In fact, if the properties of the 
Lagrange interpolation polynomials are used first, we can easily see that 



(4.26) 


/ 9 \m/2 N / 

H) TT«-v 

r , r . = s=m 

£m iy T \m/2 N / 

i 1 - 4) n <h - c s> 

s=m 


5 


where the prime indicates s ^ Jl. To transform the numerator in terms of the asso- 
ciated Legendre functions, we use the orthogonality properties of these for a fixed 
order. After some algebra, we obtain for the coefficients of equation (4. 22) 


vin/2 N / 


(N+n-m)/2 


( „\m/ z in / mni-n 

*-<?) TT«*-v- E 

s=m j=n-m 


Urn , 2j+n-m, n-m 

• . U 

2j+n-m m 


(4. 27) 


where 


^2j+n-m, n-m _ 
m 


(2n + l)(n - m)! 2 
(n + m) ! 


n! 


(n-m)/2 (2n - 2t) ! (-1)* (") 2 j_t (j + n - m - t)! 

x E — 

t=0 


(n - m - 2t) ! (2j + 2n - 2t + 1) ! ! 


and 


j0m . 

Vk = 1 

Jim 


N / 


C e 


N-k-1 " s 

s=m 


Jim 


N / N 


a N-k-2 - E E S 


s^^m s 2 =m 


4m _ . . .N-m -i— r r 
l 0 - (_1) ! 1 C s • 


N 


s=m 
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We now return to the study of the functions W^(C, M* First, we observe that for 
k = 0, since P (1) = 0, m * 0, we can write 

w oo = £ > 

n=0 

with the N + 1 requirements 

W 00<V = *0r’ r = 0, 1, • • • , N . 

This uniquely defines the sampling function Wqq(Q and identifies it with Cq(Q, de- 
fined by equation (4. 21) ; i. e. , 


= C 0 (O • 


Another important relation is found by considering that 
N 2k 

P n0«>=£ £ P nO< C k> W kj«+> 
k=0 j=0 

and also 

p„o<o = £ p„o< V c k<a - 

k=0 

so that, for k > 0, 


w 00 (0 = 


2N + 1 


INI 


1 + 2 


j=l 


T )( C) 


(4.28) 


£ W kj< C > X) =C k (C) = 2NTT 

i=o 

Recalling the form of W^.(C, given by equation (4.18) and taking equation (4. 19) 
into account, we can write 


N 

1 + 2 £ 

j=i 


Tj(0 


VV 


(4. 29) 
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T S<° - £ t S (C p ) G pm (C) 

p=m 

s S'° = 2 s S<9 V'a ■ 

p=m 

so that 

w kj (C ’ 1|S E E G pm (C> [ T S< V cos mX + ^ V sln mX ] ■ < 4 - 30> 

m=0 p=m 


The problem is now reduced to obtaining the numbers 


T S<9> s S ( V 


k = 1, 2, . . . , N ; m = 0, 1, . . . , k ; 
j = 0, 1, . . . , 2k ; p = m, m+1, . . . , N 


(4.31) 


To do so, it is convenient to consider three cases: 1) £ - £, , 
2) £ = C r , r > k. and 3) £ - £ r , r < k. 


1) Let £ = £^. We must have 


w kj< c k- V = 2 2 p - m<5k) [ A ™ 003 mX kp + sin m V] = 8 pj • 

m=0 n=m 

Multiplying successively by cos aX^ , sin aX^ and summing over p (p = 0, 1, . . . , 2k), 
we find 


N £ 

„ — , — , ,nm m 

2 P nm (C k> A kj = 5TT cos mX kj 

n=m 


N _ nTY1 9 
v p (C. ) B, = 01 ~ T sin mX, . 
2_j nm'Tc k] 2k + 1 kj 

n=m 


so that 


(4. 32) 
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(4.33) 


^.m , r , _ m , 

T kj (C k>-2irfT cos mX kj 


S S (C k> = 2iTTT SInmX kj ■ 


which, in fact, verifies equation (4.14). 


2) Consider now C = C r , r > k. We must have 


£ £ P nm<y K]“ cos mX rs + S ‘ n mX rsJ ' 0 ’ * ' k+ }’ k+2 ' • • • > N 

m=0 n=m ’ ’ ‘ ' ‘ ’ 

(4.34) 

Let m/ = 0, 1, 2, . . . , k. Then, multiplying equation (4. 34) successively by cos m' X. 


D 


nm 


,nm 


J 


rs 


and sin m' \ and finally summing over s give, respectively, 
rs 


N _ riTn 
VP (C ) Af m = 0 

/ y nm' r k] 
n=m 

N _ 

V P (C ) B? = o 

/ ,, nm'V k] 

n=m 


m = 0, 1, . . . , k ; r = k+1, k+2, . . . , N . (4. 35) 


Then, 


T S (C r> = 0 

(4. 36) 

S£<e r )=0, r>k . 

It is important to note that the results (eq. (4. 35)) are true only because in this 
case 


m + m' < 2k < 2r < 2r + 1 . 


(4. 37) 
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In another situation (r <k), the result would not be the same. In that case, the 
summation over s is not trivial, since there are divisors of the type 


sin 


m + m ; 
2r + 1 17 


(4. 38) 


that may be zero if m + m' > 2r + 1. Equations (4. 36) also show that 


w k)< C r’ X) '° . 

for any X, for r > k. 


(4. 39) 


3) Finally, let £ = £ r , r < k. In this case, the above procedure does not apply. 
Equation (4. 34) is still valid, however. Let us observe that the point distribution on 
every circle is such that 


cos mXy = cos mX k; 2k _ J+1 


sin mJy = -sin my 2k . J+1 , 


(4.40) 


j 1, 2, . . . , k , 


for any k. Therefore, substituting 2r - s + 1 for s into equation (4. 34), we get 


£ E P "” <C r ) C ° S mX rs - C sln mX rs] ' 0 > < 4 ' 41) 

m=0 n=m 


and from equations (4.34) and (4.41), we find 


k N _ nTY1 

E V P (C ) A, - cos mX = 0 
/ 4 nm'^r k] rs 

m=0 n=m 


k N _ 

E v p (C ) B, nm sin mX = 0 , 
/ 4 nm'V kj rs ’ 

m=0 n=m 


(4.42) 


or 
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(4.43) 


E T S (C r> C0S X rs = 0 

m=0 


.m 


E "S'T 

m=l 


(CJ sin X rg = 0 , 


s = 0, 1, r = 0, 1, . . . , k=l . 

Also, using equation (4.40) in equation (4.32), we find, in general, 


,nm _ .nm 
A kj A k, 2k-j+l 


nm nm 
% \ 2k-j+l ’ 


that is, 


T S (C, ' T “2k-i + l<0 


s «(« = - s ?;2k- j+ i(o 


(4. 44) 


(4. 45) 


so that the number of unknown T's and S's is reduced by half (plus one, corresponding 
to j = 0) . 


Other general important relations are obtained from 
2k 

E w kJ <CA) =C k (0 ; 

J-0 

that is, 

2k 

E T £j<o= c k® 

j=0 


2k 

^T^(Q=0 , m^O 

R 

2k 

s“(Q - 0 , m t 0 . 


(4.46) 


60 



In particular, for r/k, C, (£ ) = 0, and using equation (4.45) with equation (4.46), 

k r 

we get 


+ 2 E T S (C r> = ° 

j=l 

S^(C r )=0 , r = 0, 1, . . . , k-1 ; m = 0,l,...,k . 


(4 . 47 ) 


The system to be solved is equation (4. 43) plus equation (4. 47). But for a given pair 
(k, j), we have the system (4. 43), which alone can be solved by recurrence. In fact, 
we find 


E T ki<V + E T W <C ’ J ' ° 

m=r m=r+l 


E T U (C r' + E T k* (S r> =0 

m=r-l m=r+2 


(4. 48) 


E T W (C r ) = ° 

m=0 

for r = 0, 1, . . . , k=l, and the summations are over m modulus 2r + 1 . 

The analogous system for the S's is given by 

E S S (C r> - E ^ (C r> = 0 

m=r m=r+l 

(4. 49) 

E S H ,C r> - E ^ - 0 

m=l m=2r 

for r = 0, 1, . . . , k=l, and again the summations are over m modulus 2r + 1. 
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These two systems are obtained by recognizing that only r + 1 different values of 

COS mX can occur — that is, for example, 
sin rs ’ e ’ 

cos pX = cos (p + 2r + 1) X = . . . = cos [p + j(2r + 1)] X . 
rs rs rs 

Similar relations exist for the sines, with only an alternate change of sign. 

Finally, the values of T™(^) and needed to complete equation (4.30) are 

computed in sequence (m = k, m = k-1, . . . ), using equations (4. 33) and (4. 36) with 
(4.49). 

4.3 Approximation Error in Sampling- Function Representation 

In previous sections, we have shown that quantities generally represented by 
series of spherical harmonics can instead be expressed in a series of sampling func- 
tions. One fact, though, should always be kept in mind: Any function defined on a 
sphere enjoys the important property of least-squares deviation (and therefore con- 
vergence in the mean) when expanded in a series of spherical harmonics or in any 
other orthonormal complete set in L 2 (sphere). The coefficients of this series are 
uniquely defined. 

For truncated series, the same can be said for any sequence spanning the same 
linear manifold of the spherical harmonics. This is a basic motivation for using the 
sampling functions we have defined. 

Consider first the now familiar form of sampling functions for a sphere: 


M _ — 

2 A «w ’ 

i=l 

and, reciprocally, 

_ M _ 

x k (C ’ X) = X x k w je(CA) • 
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For spherical harmonics normalized according to equation (2. 52), we have 


_ M _ _ M _ M ___ 

= 6 kj =2 x k«i>V w r x j) W*> E A ii< x i> x j' 

«=1 4=1 i=l 


and therefore 
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X k^£>\e* A H] ’ 


M 


2^ A 4j X k ( ^’ X ^ 6 kj - 
A=1 


(4. 50) 


These relations show that the set of constants A^. is orthogonal to the set of spherical 
harmonics over the ensemble of sampling points. This property is very important, 
for if the sampling-point distribution in the limit N — 00 is everywhere dense on the 
sphere, we can conclude that 


lim 
N — * oo 


M 


E A *i x k ( W = < x j’ x k> 

i=i 




(4. 51) 


where, of course, goes into a generic point on the sphere. 

Next, consider the harmonic Fourier series for continuous functions F defined on 
the sphere: 


F(C,X) E n X n (C ’ X) ’ (4,52) 

n=l 


where = (F,X n ). The Parseval equality holds in this case; i.e. , 



On the other hand, consider the truncated representation 
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X D n (N) X n (C,X) 

n=l 


(4.54) 


From equations (4.52) and (4.54), we obtain 
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F ' f = £[ E n- D n (N) ] V C ’ X > + Y. E » X n (CA> 

n=l n=M+l 

and therefore 
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n=l n=l 
or, using the proper definitions, 
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F - f|| 2 =|| f|| 2 - 2 ^ (F, X n ) £ A kn F <Vk> 

n=l k=l 
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k=i je=i 

which, in view of equation (4. 51), leads to 
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2 £ A kn F( ^k’ X k ) A j&n F ^£ ,X ^ 


hm || F - f| 2 = I F || 2 - 2 
N — ► oo 


^(F,X n )(F,X n )+^(F,X n )(F,X n )=0 . 


n=l 


n=l 


64 



It follows that if in the limit N — ® the sampling-point distribution is everywhere 
dense on the sphere, the limit of the distance ]| F - f|| as N -* ® is zero, which 

A 

guarantees convergence in the mean of F — F. The truncation error can be obtained 
after precise bounds are found for the sampling functions or, equivalently, for the 
coefficients of their linear representation in terms of spherical harmonics. These 
have well-known estimates ( see, e.g. , Sansone, 1959, p. 202). 
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5. SAMPLING FUNCTIONS OF THREE ARGUMENTS 


5. 1 Spherical Reference Surface 

For the 1-3-5 distribution of points on a sphere or of directions in space, Section 4 
developed a set of sampling functions W (£,, X) that are linear combinations of the sur- 
face spherical harmonics. For some applications, these are the appropriate functions; 
but for others, a corresponding set of functions is needed, a set that incorporates depen- 
dence on radial distance r as well as on £ and X. A typical example is the representa- 
tion of the gravitational potential about a body such as the earth. 

Generalization to three dimensions is simple if a sphere is chosen as the reference 
figure on which the sampling points are distributed according to equations (3.3). The 
coefficients A^ in equation (4. 12) are still applicable, and the three-dimensional 
sampling functions are defined by a formula analogous to equation (4. 1), namely, 

N ^ n+1 

[ A r S nm +i rAj ’ < 5 - 1 > 

ufO m=0 

where a is the radius of the reference sphere. Conversely, 
n+l _ “ 

(f) • < 5 ' 2 > 

k=l 


As an example, the representation of the geopotential illustrates that definition 
(5. 1) achieves the desired end. The potential exterior to a body is conventionally 
written as 
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V = 
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z 

n=0 


Z I 0 ”” © n+1 X mn<^> + sMn (“f 1 

m=0 


(5.3) 


Truncating this expression at degree N, and using equation (2. 58), we can rewrite 
equation (5.3) in the form 
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k= 1 


N 


(") n+1 ^ 

Vr/ n 


ifO 
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where 
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L k = V Tc nm X K V ,U*S“Y 
n / j L nm^k’ k' irm 

m=0 


Kk’V] 


(5.4) 


(5.5) 


On the other hand, if there are (N+ l) 2 functions Z^(r, t,, X), k = 1, 2, . . . , M=(N+1) 2 , 
defined by equation (5. 1) such that, in general, 


z k<V^V = 6 k* < 5 - 6 > 

for points (r^ , , X^ ) on the sphere r^ = a, then we can write (5. 4) in the form 

M 

» < 5 - 7 > 

k=l 


where 
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k’ ^k’ V 


(5.8) 
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The form (5. 4) can be used with the right side of equation (5. 7), yielding 


*-±<sr 1 *-±& l ±F m \A, 


n=0 


r 

n=0 m=0 


X. ) + S nm Y (£, , X, ) 
k' rnn'^k’ k_ 


(5.9) 


If the spherical-harmonic coefficients C 1 ™, S 11111 are given, the corresponding 

sampling-function coefficients result from equation (5. 9). If the sampling-function 
k 

coefficients L are given, then substitution of equation (5. 1) into equation (5.7) yields 
the expressions for the spherical-harmonic coefficients: 


c nm 
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nm 
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nm 
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(5. 10) 


5. 2 Nonspherical Reference Surface 

The surface of the earth differs from a sphere by about 1 part in 298, and for some 
purposes, a sampling-point distribution on a reference ellipsoid might be preferable to 
a reference sphere. In other circumstances, a more detailed reference surface that 
better approximates an equipotenlial surface of the earth might be useful. 

* 

The objective, of course, is to formulate a set of sampling functions Z^(r, £,, X) 
defined for sampling points on the reference surface and spanning the linear manifold 
determined by the solid spherical harmonics truncated at some degree N. In principle, 
the procedures described in Section 2. 5 are applicable. In this approach, the solid 
spherical harmonics are the starting point, rather than the surface harmonics as in 
Section 4. 1. This approach maybe more or less practical, depending on how complex 
the reference surface may be and on how the points are distributed on it. 
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Another approach utilizes the Z, (r, £,X) for a spherical surface known from 

* K 

Section 5. 1. Again, the Z^ for the nonspherical reference surface must span the same 
linear manifold as do the Z^, and there must be the same number of sampling points. 

For all the applications we might immediately envision, the nonspherical figure 
will not differ vastly from a sphere, and a natural one-to-one correspondence can be 
established between sampling points on the nonspherical surface and the 1-3-5 distri- 
bution on a nearby sphere. The corresponding points can be denoted by the same 
index: 


^- k ) * * X k^ ’ k 1, 2, . . . , M (N+l) 

on a sphere on a nonspherical surface 


* 

By a sequence of transformations, the set Z^ can be carried into the set Z^ as 

follows. Consider first a set of sampling points consisting of the 1-3-5 set on a sphere, 

* (11 

except that point P^ is replaced by P^. Let Z^ ’ be the sampling functions that satisfy 
the usual conditions for this modified set of points. Then, 


M 


z k (r,S,x> = ^2 z k( r i 1) ’^ 1) ’ x i 1) ) z i 1} < r >^ x ) , 

i=l 


(5.11) 


where the sum is over the points in the modified set, associated with the zj^. But 




i + 1 • 


(5. 12) 


Thus, 


Z k (r, C, \) = (1 " 6 kl ) z[ 1} (r, £, X) + Z k (r*, £*, X*) Z™ (r, £, X) . (5. 13) 


70 



The inverse of equation (5. 13) is 


-(1) 


Zj '(r,£,X) = (1-6^) Z.(r,£,X) 

/ r * * * i * * * \ 

\ 6 il| 1 * z j^ r 1’ 1 ~ z j< r i’^i’ ) ‘i ) i 


( 


* * * 


) 


Z^^X) 


(5. 14) 


The coefficients of the Z (r, £, X) in this equation can be evaluated from the known 

1 (11 

expressions for Z^(r, £,,X). Thus, the Zj ; (r, £,,X) can now be treated as known functions 
and the same steps repeated, with a second point from the sphere being moved to the 
nonspherical reference surface. 

If the zj* 1 **(r, £,, X) are the sampling functions after h - 1 such steps, then for the 
hth step the formulas are 


z£ _1) ( r ,&,X) - (1-6^) z[ h >(r,£,X) + z[ h_1) (r*,^,X*) zf>(r,C,X) (5.15) 


and 


zj*(r, X) = (l-6 jh ) zj h_1) (r,^,X) 

6 

t 

(5. 16) 

In matrix form, the transformation (5. 16) is simply a unit matrix with the h column 
replaced by the coefficients of the second term on the right side. The successive 
application of M such transformations carries the functions Z^(r, £, X) into the functions 

z£(r,c,x)- 


jh - 


1 + z( h-1) 

_J 


* * 
r h’ ^h> 


a- 


Z (h_1) (r* t X 1 
n (r h 5 *h> n 


(h-11 * * *\ 

Z i _(h-l) 


(r, ;,X) 


71 



If the nonspherical reference surface were simple enough, the M transformations 
might be compounded analytically into a single transformation. A substantial simplifi- 
cation can also result from choosing the sampling points on the nonspherical surface 
* * 
such that ^ an d 

Of course, in choosing the new sampling points, attention must be given to 
avoiding singular cases of the transformations (5. 15) and (5. 16). Clearly, one neces- 
sary condition is that 

(h- 11 * * * 

z h ^h’V* 0 • 


For points P^ near Py this condition should be easily satisfied, since 

z h" 1 ) <v^V = 1 • 

* 

Similarly, for P^ near P^, the numerator of the second term of equation (5. 16) 
can be expected to be small if j 4 . h, since 

z r i) < r h>^V = ° • 

* 

Thus, for the case P^ near P^, transformation (5. 16) differs only slightly from an 
identity matrix. 
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6. APPLICATIONS 


6. 1 General Considerations for Applications 

As noted in Section 1, this formulation of sampling functions started in anticipation 
of the need to handle satellite-to-ocean altitude measurements. This application has 
been discussed elsewhere (Lundquist et al . , 1969; Lundquist and Giacaglia, 1972b) and 
is addressed again later in this section. It is a specific example of a larger class of 
possible geophysical applications. 

These applications occur when it is appropriate to aggregate measurements of 
some geophysical quantity into a set of representative values for localized geographical 
regions and yet maintain the ability to relate these representative values conveniently 
to a global representation of the physical quantity. For satellite-to-ocean altitudes, 
the quantity to be so treated is the geocentric radius to mean sea level. 

The sampling-function representation is, of course, formulated such that the 
tabulated values at the sampling points have a clear relation to a spherical -harmonic 
global representation. In a function-fitting sense, the dependence of many geophysical 
quantities on position can be approximated well by a spherical-harmonic expansion, 
and hence equivalently by a sampling-function expansion. Further, the spherical- 
harmonic expansion in three dimensions is manifestly appropriate for fields such as 
the gravity potential of the earth, because the spherical harmonics are solutions of 
the Laplace equation. Since the sampling functions were required to be linear com- 
binations of the spherical harmonics, they, too, are solutions of the Laplace equation. 
If this requirement had not been imposed, difficulties could arise when the gravity 
potential is expanded in functions that are not solutions of this equation. 

In the following, each specific application falls naturally into one of two groups, 
depending on whether a two- or a three-dimensional model is required. Sections 6. 2 
and 6. 3 are concerned with examples of the former group, and Section 6. 4, with the 
latter. 
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For whatever geophysical quantity is to be represented in a sampling-function 
expansion, the question naturally arises as to how values of this quantity will be cal- 
culated at points other than the sampling points. An important objective of the 
sampling-function expansion is, of course, to avoid calculating the sum of very many 
terms of a series. 


For the two-variable case of a quantity defined on a sphere or similar reference 
surface, the simplest procedure for intermediate points is the use of some elementary 
interpolation formula between the grid of values at the sampling points. For some 
quantities, there may be a particular interpolation formula that is conventionally used 
for finding values between tabulated values at a grid of points. If this is the case, the 
conventional procedure would probably be preferable. 


The general behavior of the sampling functions suggests an interpolation formula 
that can be applied if there is no reason to prefer some other particular formula. Let 
a. be the angular distance from the intermediate point to the ith sampling point. Then 
the value of a function f(£,, X) can be approximated by 


f(C,M 


E 


f(C.,X.) q 0 (a.) 


i=M 


q o< a i> 


( 6 . 1 ) 


where 


= yi - IT COS(X.-X) 


q Q (a) 


1 sin { [(2N + 1) / 2] } a 
2N+ 1 sin [( 1/ 2)a] 


from equation (2. 31), and {£2} is the set of sampling points near point (£,, X). 

The situation in three dimensions is more involved and should be considered in 
detail for specific applications. 
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6. 2 The Geoid and Mean Sea Level 


The geoid is that equipotential surface for the earth (including in the geopotential 
a term yielding centrifugal force) that best approximates mean sea level. Functionally, 
the geoid is specified when the radial distance to it from the center of mass is given as a 
function of direction — that is, when an expression for r(£, \) is provided. Thus, the 
representation of the geoid is an example in which the two-variable formulation of 
sampling functions (Section 4) can be applied. 

A convenient representation of the geoid is central to the treatment of satellite-to- 
ocean altitude measurements, because after appropriate processing, these measurements 
lead to determinations of the geocentric radius to mean sea level. Thus, a satellite 
altimeter over the open ocean traces essentially the profile of an equipotential surface 
to the extent that sea level approximates such a surface. From oceanographic studies, 
it seems that mean sea level and an equipotential surface differ by at most a few meters. 
If, however, an altimeter has an accuracy in the decimeter range, it can be expected 
to provide a representation of mean sea level to be compared with a geoid from some 
other source. 

Over land masses, information about the geoid is obtained from level measurements 
by surveying techniques and gravity measurements. This information is often compiled 
into geoid contour maps for the continental areas measured. These maps provide a 
source of data, complementary to ocean-altitude measurements, for inclusion in a 
comprehensive geoid representation. 

By means of fundamental mathematical theorems, complete knowledge of an 
equipotential surface enclosing the earth is equivalent to complete specification of the 
earth's external gravitational field, and conversely. Thus, the geopotential represen- 
tation obtained from analysis of satellite perturbations also provides information about 
the geoid — particularly about its longer wavelength features. 

In its simplest outline, a procedure for refinement of the geoid representation 
could progress as follows: First, the best available spherical-harmonic description 
of the geopotential, presumably from satellite orbit analyses, would be used to 
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produce a global geoid that could be expressed as an expansion of r(£,, X) in surface 
spherical harmonics. Second, for this geoid, the value of r(£., X.) at the sampling 
points for some suitably high value of N would be evaluated. It would also be possible 
to calculate these values directly from the three-dimensional geopotential representa- 
tion if this were computationally more efficient than the two-step procedure. Third, 
either geoid-map data for continental areas or altimeter data for ocean areas would be 
averaged or aggregated appropriately to generate refined geoid radii at the sampling 
points. Where no new data are available, the radii derived from the geopotential 
representation would be retained. Fourth, the refined geoid representation in terms 
of sampling functions could be transformed into the equivalent representation in 
spherical harmonics. Further, the geoid representation could be transformed into 
the equivalent geopotential representation. The analytical tools for the use of sampl- 
ing functions for these four steps were developed in previous sections. 

The first and fourth steps involve the auxiliary problem of the relation between 

the coefficients for an expansion of the geopotential in solid spherical harmonics and 

those for an expansion of the geoid in surface spherical harmonics. This question 

has been reviewed by Rapp (1970) and described in some detail by Bursa (1968, 1969). 

v -9 

Our discussion below follows Bursa, who excludes terms of size 10 and smaller. 


Assume that at the surface of the earth 
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The quantities a , GM, and w are assumed constant, and the normalization of the 
® * * * 
spherical harmonics X , Y is specified below. We need to adjust R^, C nm , and 

S* . 
nm 

We consider the following expressions for the radius of the geoid: 
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r = r o 22 7 P w p (C} x) 
p=i 


(6.3) 


or 
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r = R„y^ a X(£,,\) = R„> > (a X +b Y ) , 

0 / j p p VbJ ’ 0 / j / j v nm nm nm nm ’ 

p=l n=0 m=0 


(6.4) 


where 


? P - ^ r < V V 


Given r , we can compute a • and given a , we can compute r by using the formulas 
P P P P 

of Section 4. 


We note that in the two-index notation," 


r — r, . 
P kj 


a — a or b 
p nm nm 


Here, the normalization factor we adopt for X , Y (see eq. (2. 52)) is 

nnr nm - 


X = a P (C) cos m\ 
nm nm nrrr b/ 


Y = a P (£) sin m\ , 
nm nm nnr b/ ’ 
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where 


2 _ 2(2n + 1) (n - m)! 

nm £ (n + m)I 
nr ' 

(6.5) 

£ 0 = 1 ’ £ m = 2 for m > 0 » 


and we get 
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^ J dK | . 
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We note that a normalization often used (see, e.g. , Lundquist and Veis, 1966, vol. 1) 
is , 
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( 6 . 6 ) 


so that 


2tt +1 

k | dx j 

0 -1 



m= 0 
m > 0 , 


an asymmetry we should like to avoid. 


We consider first the case in which a set of coefficients a , b is given for a 

nm’ nm 

geoid representation. As a first step, we express this set in an unnormalized form, 
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a = a a 
nm nm nm 


b = a b 
nm nm nm 


(6.7) 


To calculate C , S in terms of the a , b , it is convenient to define the 
nm’ nm nm’ nm’ 

following quantities: 
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Then the parameters in a geopotential representation are 
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(eq. (6. 9) continued on next page) 
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If desired, these can be expressed in the conventional normalization: 
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(6.9) 


For the case in which the geopotential coefficients are given and the geoid coeffi- 

cients mnst be found, we first express C , S in unnormalized form: 
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The coefficients in the geoid representation are thus given by 
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Finally, these can be normalized in the convention adopted previously: 
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( 6 . 11 ) 


Relations ( 6 . 9) and ( 6 . 11) were obtained by Bursa (1969) and do not include terms 10 
or smaller. More precise relations have been derived and will be published in a 
separate article. 
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6. 3 Topography of the Earth 


The topography of the earth is another geophysical quantity for which a mathematical 
representation is sometimes needed (Lee and Kaula, 1967). This is a further case in 
which a representation in sampling functions of two arguments can be used. 

Several interrelated entities might be considered. The first is the geocentric 
radius to the mean surface of the solid earth. A second is the elevation of the mean 
solid surface above or below mean sea level. These two entities are related, of course, 
through the mean-sea-level representation discussed in Section 6. 2. 

If topographic maps of land areas and bathymetric maps of the sea floor are the 
initial sources of data, the mean elevation of the area around each sampling point 
could be averaged to give the coefficients in a sampling -function representation of the 
profile relative to mean sea level. These coefficients could then be added to the mean- 
sea-level coefficients to get the coefficients for the geocentric radius to the mean 
surface of the solid earth. Either of these representations of topography could easily 
be transformed into an equivalent spherical-harmonic representation by the formalisms 
in the previous chapters. 

The profile of the solid and ocean surface of the earth might also be desired, for 
example, as the lower boundary of the atmosphere (see, e.g., Verniani, 1966). In 
this case, the coefficients for land areas can be taken from the representation of the 
solid earth; those for the ocean areas, from the representation of mean sea level; and 
those for coastal areas, from appropriate averaging of the solid-earth and sea-level 
coefficients. Here, the advantage of the sampling-function representation is evident, 
because the required coefficients can be selected individually on the basis of the geo- 
graphical situation in the neighborhood of the sampling points. 

6.4 The Geopotential 


The gravitational potential around the earth is an important example of a quantity 
that requires description in functions of three arguments. The formalism of Section 5 
was presented with this in mind. Section 6. 2 discussed the important relations between 
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the geoid and the geopotential. If a spherical reference surface of radius R is used, 
following the formalism of Section 5. 1, the potential can be written 
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or equivalently as 
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Then the coefficients are related by 
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and inversely, 
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(6. 15) 


with the and B^ m as derived in Section 4. Thus, a sampling-function repre- 
sentation can be used effectively whenever it may be a convenient alternative to an 
expansion in spherical harmonics. ’ 

Other alternatives to a spherical-harmonic expansion have been advanced in the 
geophysical literature, and their relationship to a sampling-function representation 
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can be examined. Several methods employ adjustable mass distributions on a spherical 
shell to achieve a detailed representation of the geopotential (e.g. , Vinti, 1971; 
Croopnick, 1970; Koch and Morrison, 1970; Koch and Witte, 1971; Koch, 1972). Nor- 
mally, the leading inverse-square term, the oblateness term, and perhaps a few other 
low-degree terms would be represented in the usual spherical -harmonic form, and the 
surface mass distribution would represent only the remaining short-wavelength detail 
in the potential. 

The basic theoretical idea behind this method is the consideration that a spherical 
shell of surface mass density 
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produces an external potential of the form (6. 12), where 
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(6. 17) 


and as usual. 
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(6. 18) 


and similarly for Y . 

nm 

If the spherical surface is divided into a grid of areas S^, the surface integral 
(6. 18) can be written as a sum of integrals over these areas. For simplicity, the 
density in each area can be assumed to be a constant, Then, 
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(6. 19) 


4it(x ’ X nm ) x k X nm (S k ) ’ 

k 

where 

wv = I ^ dS ■ < 6 - 20 > 

s k 

The integrals X ^ m (S^) depend only on the spherical harmonics and on the way the 
sphere is divided into the areas S^. These integrals are constants that can be easily 
calculated and tabulated. 

Substituting equation (6. 19) into equation (6. 17) gives the coefficients for the 
spherical-harmonic expansion of the potential due to the specified mass distribution x k ' 
This will be an infinite set of coefficients, which, of course, can be suitably truncated. 
However, a desired objective is to avoid calculating the spherical harmonics; this can 
be accomplished by an approximation that replaces each area of surface mass by a 
centered point mass (Sjogren et al . , 1971; Balmino, 1972; Kaplan and Kunciw, 1972; 
Needham, 1970), 


m k = *k S k 


The potential of such a distribution of point masses is 


V= - 



> 


( 6 . 21 ) 


( 6 . 22 ) 


where 


= (x -x k ) 2 + (y-y k ) 2 + (z - z k ) 2 , 

( x k’ x k’ z k^ being the coordinates of the point masses. 
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If one should wish to adopt this sort of formalism, use of the sampling functions 
can still effect some simplifications. The first step is to recognize that the wavelength 
of the meaningful detail in the potential is limited by the degree at which equation (6. 16) 
is truncated, by the spacing of the areas S^, or by the spacing of the point masses. 
Thus, it is appropriate to pick the degree N corresponding to the desired short -wave- 
length limit and to work within the linear m anif old defined by the spherical harmonics 
through degree N. 

The second step is to expand x(C,X) i n sampling functions: 

M 

= • < 6 - 23 ) 
i=l 


Given a set of coefficients C 11111 , S 11111 , the values of x (?>p X.) follow simply from the 
evaluation of equation (6. 16) at the sampling points. Presumably, the values x(C i , X.) 
would be subject to an adjustment that makes use of some sort of observational data. 

The next step is to find the coefficients C nm , S 1 ™ that correspond to the improved 
X (^, X.). Substituting equation (6. 23) into equation (6. 16) gives 

2 M 

cmn - T, X «i> X i> ( W i«> X), X nm K, X)) . <6. 24) 

i=l 

But 

(Wift, X), X^, X)) = A* m , (6. 25) 
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so 


-nm 


2 M 

4ttR \ A . , A nm 
(2n+ 1) p 2 If x ^i’V A i 
i=l 


s nm 


4ttR 2 
(2n+ 1) |i 


M 


i=l 


’ r l 


nm 


(6.26) 


These are the desired equations for calculating the potential coefficients from the 
model values of the surface density. 

In summary, if it is desired to use a surface -density model for the geopotential, 

2 

it is necessary to use only the (N+ 1) density values at the defined sampling points as 
the model parameters. The geopotential coefficients for a spherical-harmonic expan- 
sion (or a sampling-function expansion) follow directly from those values. 

It is also possible to take the further step to a mass-point representation. The 
area associated with each sampling point (t,^, X.^) can easily be specified. A mass 
m^ can hence be associated with the point by 

m k = f X(&,MdS , (6.27) 

S k 

or, with lesser accuracy, by 




m k x( ^k’V S k ' 


(6.28) 


These can be used in a potential expression of the form of equation (6. 22). 
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6. 5 Other Applications 


Many other applications of sampling-function formalism could be discussed. 
Perhaps the representation of vector fields is the most important. Nevertheless, 
these other applications, in general, parallel the examples of the previous sections 
and thus will not be treated here. We leave such discussions until they are needed 
in conjunction with a specific application. 
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